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1.  Introduction 

v'  v  r> 

This  project  deals  with  the  problem  of  solving  large  systems  of  stiff  ODF«  In  particular,  we  see  as  one  of  the 

major  issues.the  choice  of  methods  for  solving  the  systems  of  linear  and  nonlinear  equations  that  arise  at  each 
I  *  1  s  f 

integration  step.  -Wr  proposed  to  use  variants  of  the  preconditioned  conjugate  gradient  (PCG)  method  for 

solving  the  linear  equations,  and  truncated  Newton-like  iterations  for  solving  the  nonlinear  equations.'  Our  main 

objective  is  to  carry  out  a  systematic  study  of  the  methods  involved  and  ultimately  to  produce  a  well- 

documented  and  well-tested  computer  code  for  solving  large  systems  of  stiff  ODBs  that  incorporates  the  results 

of  our  study.  A  brief  summary  of  our  original  proposed  research  will  be  given  in  Section  2. 

In  the  17  months  (1  year  plus  5  months  no-cost  extension)  supported  by  this  grant,  we  believe  we  have  made 
definite  progress  towards  our  ultimate  goal.  We  have  developed  some  fundamental  algorithms  in  the  area  and 
produced  a  first  version  of  the  program  package  that  incorporates  PCGPACK  with  LSODE.  Highlights  of  the 
results  accomplished  will  be  presented  in  Section  3.  A  total  of  five  technical  reports  have  been  produced  under 
this  grant,  of  which  three  have'been  accepted  for  publications  in  SIAM  journals,  and  one  is  near  completion. 
Copies  of  thes^  reports  have  already  been  submitted  to  AFOSR  before  and  will  not  be  included  with  this  final 
report. 

2.  Goals  of  the  Project 

Our  main  objective  is  to  determine  whether  the  use  of  iterative  methods  for  solving  linear  systems  of  algebraic 

equations  in  codes  for  large  stiff  systems  of  ODEs  is  competitive  with  sparse  direct  techniques.  A  more  tangible 

objective  is  to  produce  a  computer  program  that  incorporates  such  an  approach.  There  are  at  least  four  main 

issues  that  have  to  be  investigated: 

.  & 

1.  Which  preconditioned  conjugate  gradient  method  is  the  best  for  our  purpose?  (Which  variants  of 
the  method:  Conjugate  Gradient  or  Conjugate  Residual?  Which  of  the  many  nonsymmctric 
variants?  Which  preconditionings?) 

2.  What  is  the  best  way  to  solve  the  nonlinear  systems?  (Newton-like  methods,  truncated  versions, 
chord  versions,  directional  differencing,  efficient  evaluation  of  the  Jacobian.) 

3.  How  to  modify  the  formula,  strategies,  and  heuristics  of  the  ODE  solver  to  take  advantage  of  the 
greater  freedom  offered  by  the  iterative  methods?  (Variable  coefficient  formulas,  more  frequent 
change  of  step  size  and  order.) 

4.  How  does  a  code  based  on  our  approach  compare  with  the  existing  codes  in  its  performance  on  real 
problems?  (Electrical  network  problems,  Method  of  Lines  for  PDEs.) 


3.  Highlights  of  Results 


3.1.  Alternating  Direction  Incomplete  Factorisations 

The  use  of  a  good  preconditioning  is  often  essential  for  the  successful  application  of  the  conjugate  gradient 
method.  One  of  the  most  promising  class  of  preconditionings  is  the  so-called  Incomplete  Factorizations, 
which  can  be  viewed  as  factoring  the  coefficient  matrix  approximately  under  the  constraint  that  the  factors  have 
to  be  sparse.  For  elliptic  PDEs,  quite  a  lot  of  theory  has  been  developed  for  these  preconditionings.  Most  of 
these  preconditionings  have  the  property  that,  when  applied  to  the  discretization  of  a  self-adjoint  elliptic 
opertor,  the  number  of  iterations  required  to  reduce  some  norm  of  the  error  by  a  factor  of  t  is  0(h'*/2log(l/()). 
where  h  is  the  mesh  size  used  in  the  discretization.  By  modifying  one  of  these  preconditionings,  namely,  the 
Dupont-Kendall-Rachford  preconditioning  (DKR),  we  have  been  able  to  derive  some  error  estimates  which 
suggest  that  the  number  of  iterations  can  be  reduced  to  0(h‘1/,3log(l/c)).  Thus  this  new  preconditioning,  which 
we  call  ADDKR,  is  asymptotically  faster  than  DKR.  Moreover,  numerical  experiments  on  some  model  elliptic 
problems  indicate  that  even  for  problems  of  moderate  size  the  new  preconditioning  is  competitive  with  DKR. 
We  feel  that  this  is  an  important  discovery,  both  theoretically  and  practically,  since  no  other  known 
preconditionings  in  this  class  have  been  shown  to  possess  a  better  asymptotic  rate  of  convergence  than  DKR. 
This  work  appeared  in  the  SIAM  Journal  of  Numerical  Analysis  in  April,  1983. 

3.2.  Nonlinear  Preconditionings 

One  promising  idea  that  some  people  have  noticed  in  the  application  of  Krylov  subspace  methods  for  solving 
the  linear  systems  Aw  =  b  that  arise  in  the  inner  loop  of  a  Newton-like  iterative  method  for  solving  nonlinear 
systems  F(x)  —  0  (with  A  being  the  Jacobian  of  F  at  a  point  x)  is  the  use  of  the  directional  differencing  ( 
F(x  +  cu)  *  F(x)  )  /  c  for  approximating  the  matrix-vector  product  Au  in  the  CG  algorithm.  This  avoids  explicit 
evaluations  of  the  Jacobian  and  only  requires  function  evaluations.  However,  most  of  the  preconditioning 
techniques  in  use  now  are  derived  from  the  matrix  elements  explicitly.  It  is  therefore  unclear  as  to  how  to  apply 
the  preconditionings  when  the  matrix  is  not  explicitly  available.  This  issue  does  not  seem  to  have  been 
addressed  in  the  literature.  We  have  obtained  some  results  in  this  direction.  We  have  derived  an  algorithm  for 
preconditioning  the  CG  iteration  with  directional  differencings  that  reduces  to  the  SSOR  preconditioning  in  the 
linear  case.  It  only  requires  evaluating  the  diagonal  elements  of  the  Jacobian  which  can  also  be  approximated 
by  function  evaluations  and  can  be  easily  incorporated  into  the  CG  iteration.  We  consider  it  an  important 
discovery  because  we  think  the  use  of  directional  differencing  will  prove  to  be  very  useful  in  the  stiff  ODE 
context.  We  plan  to  incorporate  the  results  into  the  final  version  of  our  code.  We  presented  a  talk  on  this  topic 
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in  the  Sparse  Matrix  Symposium  held  in  Fairfield  Glade,  Tennessee  in  October,  1982.  The  paper  has  been 
accepted  for  publication  in  the  SIAM  Journal  of  Scientific  and  Statistical  Computing. 

3.3.  Basie  CG  Variant  of  LSODE 

All  versions  of  our  code  are  going  to  be  built  around  and  incorporated  into  the  package  LSODE  by  Alan 
Hindmarsh.  A  first  version  of  the  code  has  been  completed.  It  basically  replaces  the  sparse  matrix  solver  in 
LSODE  with  the  conjugate  gradient  solver  PCGPACK  developed  at  Yale.  The  code,  with  we  call  LSODECG, 
can  now  take  problems  in  the  sparse  format  (same  as  the  IA-JA  format  used  in  the  Yale  Sparse  Matrix  Package 
(YSMP)  and  in  LSODES)  and  solves  the  linear  systems  that  arise  at  each  integration  step  by  variants  of  the 
preconditioned  conjugate  gradient  method.  It  computes  the  Jacobian  exactly  by  rows  and  uses  absolute  error 
control.  One  of  the  more  subtle  issues  is  the  control  of  the  various  levels  of  iterations  in  the  code:  the  ODE 
step,  the  nonlinear  iteration  at  each  step  and  the  linear  CG  iteration  at  each  nonlinear  iteration.  For  example, 
the  stopping  criterion  for  the  conjugate  gradient  method  has  to  be  related  to  the  user  supplied  error  tolerance 
for  the  ODEs.  We  have  obtained  some  very  useful  results  in  this  direction. 

Another  unexpected  discovery  is  that  in  solving  the  linear  system  (I  -  h/?J)x  =  b  at  each  integration  step, 
where  x  is  the  change  in  the  Newton  iterates,  it  seems  to  be  better  to  take  as  initial  guess  for  x  the  right  hand 
side  b  rather  than  the  natural  choice  of  the  zero  vector.  An  intuitive  explanation  for  this  is  that  the  coefficient 
matrix  is  approximately  equal  to  the  identity  matrix  in  the  subspace  of  the  small  eigenvalues  of  J  whereas  the 
component  of  x  in  the  dominant  subspace  is  small  outside  the  transient  regions.  Therefore,  the  right  hand  side 
b  is  a  good  approximation  to  the  solution  x.  This  has  very  important  practical  implications  since  it  appears  to 
reduce  the  number  of  CG  iterations  greatly. 

We  have  performed  some  preliminary  tests  on  LSODECG  with  problems  in  the  STIFF  DETEST  program 
(developed  at  the  University  of  Toronto  for  testing  stiff  ODE  solvers)  and  on  some  model  PDEs.  Our  main 
objective  is  to  compare  it  the  sparse  option  LSODES  in  the  LSODE  package.  These  preliminary  results  seem  to 
confirm  the  effectiveness  of  the  use  of  CG-like  methods,  especially  for  larger  problems  like  those  arising  from 
3-D  PDEs.  For  simple  2-D  PDEs,  however,  it  seems  difficult  to  beat  LSODES.  This  is  because,  for  smooth 
problems  where  LSODES  can  get  away  with  frequently  not  refactoring  the  Jacobian  matrix  and  reusing  an  old 
factorization  in  a  chord  Newton  iteration,  the  factorization  cost  becomes  negligible  and  this  makes  LSODES 
very  efficient  on  such  problems.  We  expect  the  situation  will  be  more  in  LSODECG's  favor  for  more  difficult 

problems  where  the  Jacobian  matrix  changes  more  frequently. 

I 

A  report  on  these  results  is  under  preparation  and  should  be  available  towards  the  end  of  summer,  1983. 


3.4.  Lanczos  Method  for  Multiple  Right  Hand  Sides 

In  the  stiff  ODE  context,  it  is  important  to  be  able  to  deal  with  systems  with  multiple  right  hand  sides.  This 
is  because  the  right  hand  sides  are  usually  not  available  simultianeously  as  they  depend  on  the  solution  of  the 
previc**s  system.  One  way  of  achieving  improved  efficiency  is  by  saving  the  Krylov  subspaces  generated  for  the 
previous  systems  and  reusing  them  in  subsequent  solutions,  as  proposed  by  Parlett.  We  have  obtained  some 
results  in  this  direction.  We  showed  that  one  of  the  algorithms  proposed  by  Parlett  is  theoretically  equivalent  to 
the  Block  Lanczos  method.  A  technical  report  on  this  topic  has  been  completed. 

3.5.  Comparison  of  Krylov  subspace  methods 

A  host  of  algorithms  belonging  to  the  same  class  of  Krylov  subspace  methods  have  been  recently  at  the  focus 
of  several  researchers  dealing  with  the  solution  of  large  sparse  unsym metric  systems.  Among  them  we  mention 
the  subclass  of  the  ORTHOMIN(kXVinsome),  GCR(k)  algorithms,  (Eisenstat,  Elman  and  Schultz). 
ORTHORES(k),  ORTHODIR(k)  developed  by  Young  and  Jea,  IOM(k)  (Saad),  Axelsson's  method,  Lanczos" 
method  and  others.  Although  the  literature  is  rich  in  methods,  few  comparisons  of  the  performances  have  been 
proposed.  In  the  particular  context  of  stiff  ODE  the  choice  of  a  method  is  even  less  clear. 

In  a  recent  paper,  we  have  performed  some  comparisons  with  a  class  of  Krylov  subspace  methods,  both 
theoretically  and  empirically.  In  particular,  it  was  shown  that  the  ORTHORES(k)  method  of  Jea  and  Young  is 
in  fact  equivalent  to  the  IOM(k)  method,  i.e.  the  iterates  produced  by  both  algorithms  are  identical.  On  the 
other  hand  both  numerical  stability  and  computational  cost  are  in  favor  of  the  IOM(k)  algorithms.  A  number  of 
comparisons  shown  in  Howard  Elman's  recent  PhD  thesis  at  Yale  favor  the  Chebychev  based  algorithms,  in 
general,  against  the  Conjugate  Gradient  type  methods.  A  common  argument  raised  against  the  Chebychev 
algorithm  is  that  it  requires  the  knowledge  of  some  of  the  eigenvalues  of  A.  While  there  is  no  doubt  that  the 
process  of  estimating  the  eigenvalues  and  obtaining  the  best  ellipse  containing  the  spectrum  of  A  may  seem 
somewhat  cumbersome,  there  are  a  few  important  points  that  can  make  the  process  highly  competitive: 

1.  Once  the  optimal  parameters  are  known,  the  method  is  extremely  fast, 

2.  Hybrid  methods  combining  conjugate  gradient  type  method  to  estimate  eigenvalues  together  with 
subsequent  use  of  Chebychev  type  methods  may  turn  out  to  be  very  efficient, 

3.  In  the  stiff  ODE  context,  the  problem  of  estimating  the  parameters  is  eased  by  the  fact  that  they 
are  likely  to  evolve  smoothly  thus  allowing  the  use  of  a  set  of  parameters  for  a  large  number  of 
integration  steps. 

Although  the  choice  of  the  best  iterative  method  in  the  stiff  ODE  context  is  still  far  from  resolved,  and  it  may 
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never  be,  we  feel  that  we  have  a  much  better  understanding  of  the  relative  performance  of  the  various 
algorithms  now.  The  paper  describing  this  work  has  been  accepted  to  appear  in  the  SIAM  Journal  of  Sci.  St  at. 
Comp. 
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